#include <bits/stl_alloc.h>
#include <stdexcept>
#include <string>
#include <algorithm>
#include <iostream>
#include <vector>
#include <sstream>
#include <deque>
#include <boost/lexical_cast.hpp>
#include <boost/tuple/tuple.hpp>
#include "bl2seq_mine.h"
#include "portable_timer.h"
#include "batch_job.h"
#include "progress_indicator.h"
#include "pg_copy_guard.h"
#include "pg_db_cursor.h"
#include "get_arg.h"
# define DEBUG_LEVEL 9
#include "debug_func.h"
#include "stlhelper.h"
#include "string_hash.h"
#include <ext/hash_map>
#include <gc.h>

using std::cerr;
using std::cout;
using std::string;
using std::vector;
using std::deque;
using std::sort;
using std::ostringstream;
using boost::lexical_cast;
using boost::tie;

//________________________________________________________________________________________

//		data structures, hashes etc.

//________________________________________________________________________________________

struct t_align_error
{
	string gene_id1;
	string gene_id2;
	string error_str;
	t_align_error(const string& g1, const string& g2, const string& err)
			:gene_id1(g1), gene_id2(g2), error_str(err){}
};
	
struct t_gi_seq
{
	string gene_id;
	string acc;
	string pi;
	string seq;
	t_gi_seq(const string& g, const char* a, const char* p, const char* s):
				gene_id(g), acc(a), pi(p), seq(s)
	{
		bl_prepare_sequence_str(seq);
	}
};

inline std::ostream& operator << (std::ostream& ostrm, const t_gi_seq& prot)
{
	ostrm
			<< prot.gene_id		<< "\t"
			<< prot.acc    		<< "\t"
			<< prot.pi;
	return ostrm;
}


typedef deque<t_gi_seq> t_array_gi_seq;
// hash map
typedef __gnu_cxx::hash_multimap<string, unsigned, hash_string> t_map_gi_index;
typedef t_map_gi_index::const_iterator t_map_gi_index_iter;






//________________________________________________________________________________________

//		save_alignment_to_db

//________________________________________________________________________________________
void save_alignment_to_db(	t_pg_copy_guard&		to_db,
							const t_gi_seq&			prot1,
							const t_gi_seq&			prot2,
							t_blast_results&		results)
{
	debug_func_cerr( 8 );
	ostringstream ostrm;
	ostrm	<< prot1	<< "\t"
			<< prot2	<< "\t"
			<< results	<< "\n";
	to_db += ostrm.str();
	to_db.cached_write();
}
	

//________________________________________________________________________________________

//		get_best_aligning_sequence_of_gene

//			get best alignment among all vs all specified proteins in 2 lists
//			usually for retrieving best aligning sequence of gene

//________________________________________________________________________________________
bool get_best_aligning_sequence_of_gene(const t_array_gi_seq& 	array_gi_seq,
										std::vector<unsigned>	vec_prot1,
										std::vector<unsigned>	vec_prot2,
										t_blast_options& 		blast_options,
										t_blast_results&		best_results,
										unsigned&				best_index1,
										unsigned&				best_index2,
										unsigned& 				cnt_blast_comparisons)
{
	debug_func_cerr( 8 );
	bool aligned = false;
	// best_results.error_message;
	best_results.bitscore = 0;
	unsigned min_seq_len = 60000;
	for (unsigned i = 0; i < vec_prot1.size(); ++i)
	{
		for (unsigned j = 0; j < vec_prot2.size(); ++j)
		{
			++cnt_blast_comparisons;
			min_seq_len = std::min(min_seq_len, array_gi_seq[vec_prot1[i]].seq.length());
			min_seq_len = std::min(min_seq_len, array_gi_seq[vec_prot2[j]].seq.length());
			t_blast_results curr_results;
			curr_results.error_message = "No similarity";
			if (!do_bl2seq(	blast_options,
							array_gi_seq[vec_prot1[i]].seq,
							array_gi_seq[vec_prot2[j]].seq,
							curr_results))
			{
				continue;
			}

			if (curr_results.align_len < 50)
			{
				// save this error if nothing else has aligned
				if (!aligned)
					best_results.error_message += ":Too short";
				continue;
			}
			if (curr_results.e__value > 1e-5)
			{
				// save this error if nothing else has aligned
				if (!aligned)
					best_results.error_message += ":e_value too big";
				continue;
			}
			
			//
			//	Is this a better alignment?
			//
	
			// save only if bitscore larger
			if (curr_results.bitscore < best_results.bitscore)
				continue;

			// if same bitscore
			else if (curr_results.bitscore == best_results.bitscore)
			{
				//only save if lower e_value
				if (curr_results.e__value > best_results.e__value)
					continue;	
	
				//  or shorter length
				if (curr_results.total_len  >= best_results.total_len)
					continue;
			}
	
			//
			//	Save better alignment
			//
			best_results	= curr_results;
			best_index1		= vec_prot1[i];
			best_index2		= vec_prot2[j];
			aligned			= true;
		}
	}
	if (!aligned)
	{
		best_results.error_message  = "No similarity. "
								"Shorter seq = " +
								lexical_cast<string>(min_seq_len);
	}
	debug_cerr(8, "get_best_aligning_sequence_of_gene normal exit");
	return aligned;

}


//________________________________________________________________________________________

//		get_all_proteins_for_gene

//________________________________________________________________________________________
bool get_all_proteins_for_gene(	std::vector<unsigned>&		vec_prot,
								const t_map_gi_index&		map_gi_index,
								const std::string&			gene_id)
{
	t_map_gi_index_iter first, last;
	tie(first, last) = map_gi_index.equal_range(gene_id);
	for (t_map_gi_index_iter i = first; i != last; ++i)
		vec_prot.push_back(i->second);
	return vec_prot.size() != 0;
}	

//________________________________________________________________________________________

//		dump_align_fail_message

//________________________________________________________________________________________
void dump_align_fail_message(std::ostream& 					BLAST_LOG,
							const std::vector<unsigned>&	vec_prot1,
							const std::vector<unsigned>&	vec_prot2,
							const string&					gene_id1,
							const string&					gene_id2,
							const t_array_gi_seq& 			array_gi_seq,
							const string&					error_message)
{
	BLAST_LOG	<< gene_id1 << "\t"	
						<< gene_id2 << "\t"
						<< error_message << "\n";
	BLAST_LOG << string(80, '-') << "\n";
	for (unsigned i = 0; i < vec_prot1.size(); ++i)
		BLAST_LOG <<  bl_retrieve_string_from_prepared_sequence(
									array_gi_seq[vec_prot1[i]].seq) << "\n";
	BLAST_LOG << string(80, '-') << "\n";
	for (unsigned i = 0; i < vec_prot2.size(); ++i)
		BLAST_LOG <<  bl_retrieve_string_from_prepared_sequence(
									array_gi_seq[vec_prot2[i]].seq) << "\n";

	BLAST_LOG << string(80, '-') << "\n";
}


//________________________________________________________________________________________

//		align_gene_pair

//________________________________________________________________________________________
bool align_gene_pair(t_pg_copy_guard& 		to_db,
					const t_array_gi_seq& 	array_gi_seq,
					const t_map_gi_index& 	map_gi_index,
					string					gene_id1,
					string					gene_id2,
					t_blast_options& 		blast_options_seg,
					t_blast_options& 		blast_options,
					unsigned& 				cnt_blast_comparisons,
					unsigned& 				cnt_aligned_pairs,
					t_program_args&			args,
					deque<t_align_error>&	errors)
{
	debug_func_cerr( 8 );

	// make sure save with gene_id1 < gene_id2
	if (gene_id1 > gene_id2)
		std::swap(gene_id1, gene_id2);
	
	
	// get protein indices for all transcripts for the gene pair
	//
	std::vector<unsigned> vec_prot1;
	if (!get_all_proteins_for_gene(vec_prot1, map_gi_index, gene_id1))
	{
		errors.push_back(t_align_error(gene_id1, gene_id2, "Missing sequence for " + gene_id2));
		return false;
	}
	std::vector<unsigned> vec_prot2;
	if (!get_all_proteins_for_gene(vec_prot2, map_gi_index, gene_id2))
	{
		errors.push_back(t_align_error(gene_id1, gene_id2, "Missing sequence for " + gene_id2));
		return false;
	}
	
	
	t_blast_results		blast_results;
	unsigned best_i = 99999;
	unsigned best_j = 99999;

	
	// Find any seg-ed blast
	if(!get_best_aligning_sequence_of_gene(	array_gi_seq, vec_prot1, vec_prot2,
										blast_options_seg, blast_results,
										best_i, best_j,cnt_blast_comparisons))
	{
		// Otherwise find any blast align without seg
		if(!get_best_aligning_sequence_of_gene(	array_gi_seq, vec_prot1, vec_prot2,
										blast_options, blast_results,
										best_i, best_j,cnt_blast_comparisons))
		{
			errors.push_back(t_align_error(gene_id1, gene_id2, blast_results.error_message));
			//dump_align_fail_message(args.BLAST_LOG(), vec_prot1, vec_prot2,
			//						gene_id1, gene_id2, array_gi_seq, blast_results.error_message);
			return false;
		}
	}
	else
	{
		//
		//	Find un'seg'ed alignment of best seg-ed alignment
		//
		//	Just test these two proteins
		vec_prot1.clear();
		vec_prot1.push_back(best_i);
		vec_prot2.clear();
		vec_prot2.push_back(best_j);
		
		t_blast_results		blast_results_seg = blast_results;
		blast_results = t_blast_results();
		if(!get_best_aligning_sequence_of_gene(	array_gi_seq, vec_prot1, vec_prot2,
										blast_options, blast_results,
										best_i, best_j,cnt_blast_comparisons))
		{
			string error_str = blast_results.error_message +
							" (without seg. Alignment possible with seg!!)\t";
			//dump_align_fail_message(args.BLAST_LOG(), vec_prot1, vec_prot2,
			//						gene_id1, gene_id2, array_gi_seq, error_str);
			save_alignment_to_db(to_db, array_gi_seq[best_i],
										array_gi_seq[best_j],
										blast_results);
			return true;
		}
	}
	
	save_alignment_to_db(to_db, array_gi_seq[best_i],
								array_gi_seq[best_j],
								blast_results);
	++cnt_aligned_pairs;
	debug_cerr(8, "align_gene_pair normal exit\n");
	return true;

}


//________________________________________________________________________________________

//		align_all_pairs

//________________________________________________________________________________________
bool align_all_pairs(	t_pg_cursor& 			pg_db,
						const t_array_gi_seq& 	array_gi_seq,
						const t_map_gi_index& 	map_gi_index,
						const t_gi_pairs&		gi_pairs,
						t_blast_options& 		blast_options_seg,
						t_blast_options& 		blast_options,
						unsigned& 				cnt_blast_comparisons,
						t_program_args&			args,
						deque<t_align_error>&	errors
						)
{
	debug_func_cerr( 8 );
	t_pg_copy_guard copy_buffer(pg_db, args.output_table +
								"(gene_id1,acc_code1, pi1, gene_id2, acc_code2, pi2," +
								t_blast_results::get_sql_column_names() + ")", args.CERR());
	if (!copy_buffer)
		return false;

	args.VERBOSE() << "Aligning " << gi_pairs.size() << " genes pairs ...\n";
	unsigned cnt_aligned_pairs = 0;
	
	{
		t_progress_indicator dots(args.VERBOSE(), std::min(20U, gi_pairs.size() / 50), gi_pairs.size());
		dots.use_timer();
	
		
		for (unsigned i = 0; i < gi_pairs.size(); ++i)
		{
			align_gene_pair(	copy_buffer,
									array_gi_seq ,
									map_gi_index,
									gi_pairs[i].gene_id1,
									gi_pairs[i].gene_id2,
									blast_options_seg,
									blast_options,
									cnt_blast_comparisons,
									cnt_aligned_pairs,
									args,
									errors);
			++dots;
		}
	}
	// save results
	if (!copy_buffer.end_copy())
		return false;
	
	args.VERBOSE() << cnt_aligned_pairs << " genes pairs aligned.\n";
	return true;
}





		
		
		
		
		
//________________________________________________________________________________________

//		get_gene_sequences

//________________________________________________________________________________________
bool get_gene_sequences(t_pg_cursor& pg_db,
						const string& batch,
						t_array_gi_seq& array_gi_seq,
						t_map_gi_index& map_gi_index,
						t_program_args& args)
{
	debug_func_cerr( 8 );
	string sql_cmd = "SELECT "
							"DISTINCT ON (a.gene_id, acc_code) "
							"a.gene_id, acc_code, s.pi, sequence "
					 "FROM "
							"sequences s, "
							+ args.map_table + " e, "
							"("
								"SELECT "
									"gene_id1 AS gene_id "
								"FROM "
									"batch_gene_pairs "
								"WHERE "
									"batch = " + batch + " AND "
									"category = '" + args.job_category + "' "
								"UNION "
								"SELECT "
									"gene_id2 AS gene_id "
								"FROM "
									"batch_gene_pairs "
								"WHERE "
									"batch = " + batch + " AND "
									"category = '" + args.job_category + "' "
							") as a "
						"WHERE "
							"s.pi = e.pi AND "
							"e.gene_id = a.gene_id";
	// create binary cursor
	if (!pg_db.Declare(sql_cmd, false))
		return pg_print_sql_err("create a cursor to "
								"retrieve sequences for batch =[" + batch + "] and "
								"category = ['" + args.job_category + "']",
								pg_db,
								sql_cmd, args.CERR());
	
	{
		unsigned batch_size = 500;
		// get data in batches
		while (pg_db.Fetch(batch_size))
		{
			int curr_batch_size = pg_db.Tuples();
			if (!curr_batch_size)
				break;
			for (int i = 0; i < curr_batch_size; ++i)
			{
				// save gi and sequence etc in deque which grows efficiently
				assert(pg_db.getvalue(i, 0));
				assert(pg_db.getvalue(i, 1));
				assert(pg_db.getvalue(i, 2));
				assert(pg_db.getvalue(i, 3));
				string gene_id = pg_db.getvalue(i, 0);
				array_gi_seq.push_back(t_gi_seq(	gene_id,
													pg_db.getvalue(i, 1),
													pg_db.getvalue(i, 2),
													pg_db.getvalue(i, 3)));
				// save deque position in a hash for fast lookup
				// N.B. hash is of c_str which is OK only for deque which doesn't
				// reallocate memory. Otherwise the pointer would be invalidated
				map_gi_index.insert(t_map_gi_index::value_type(gene_id,
																array_gi_seq.size() - 1));
			 }
	
		}
	
	}
	
	// close cursor
	if (!pg_db.Close())
		return pg_print_sql_err("close the  cursor for "
								"retrieving sequences for batch =[" + batch + "] and "
								"category = ['" + args.job_category + "']",
								pg_db, "CLOSE CURSOR", args.CERR());
	return true;
}




//________________________________________________________________________________________

//		write_alignment_errors

//________________________________________________________________________________________
bool write_alignment_errors(t_pg_cursor& pg_db, deque<t_align_error>&	errors, t_program_args& args)
{
	debug_func_cerr( 5 );
	if (errors.size() == 0)
		return true;
	t_pg_copy_guard copy_buffer(pg_db, args.error_table +
								"(gene_id1, gene_id2, message)", args.CERR());
	if (!copy_buffer)
		return false;

	for (unsigned i = 0; i < errors.size(); ++i)
	{
		copy_buffer +=	errors[i].gene_id1 + "\t" +
						errors[i].gene_id2 + "\t" +
						errors[i].error_str + "\n";
	}
	// save results
	if (!copy_buffer.end_copy())
		return false;
	
	return true;
}


	




//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888






void do_nothing_for_gc_warnings(char */*msg*/, GC_word /*arg*/)
{
	debug_func_cerr(6);
//	if (strncmp("GC Warning: Repeated allocation", msg, 31) == 0)
//		return;
//	cerr << msg << "\n";
}
	


int main(int argc, char* argv[])
{
	GC_set_warn_proc(do_nothing_for_gc_warnings);
	
	debug_func_cerr( 8 );
	
	// don't use C streams
	std::ios::sync_with_stdio(false);

	// get command line arguments
	t_program_args args;
	if (!GetArguments(argc, argv, args))
		return 1;

	
	//____________________________________________________________________________________
	// 		Connect to PostgreSQL
	//____________________________________________________________________________________
	t_pg_cursor pg_db(args.db_login.c_str(), "panda_cursor");
	if ( pg_db.ConnectionBad() )
	{
		cerr	<< "\nERROR!! Connection to database ["
							<< args.db_login << "] failed.\n"
							<< pg_db.ErrorMessage() << "\n\n";
		return 1;
	}

	args.open_log_files(get_backend_pid(pg_db));


	unsigned cnt_batches = 0;
	t_portable_timer timer;
	t_blast_options blast_options_seg(true, 1e-5, "BLOSUM80");
	t_blast_options blast_options(false, 1e-5, "BLOSUM80");
	unsigned cnt_blast_comparisons = 0;
	
	pg_db.Commit();

	
	try
	{
	
	t_batch_num(pg_db).free_hanging_jobs();
	while (1)
	{
		//________________________________________________________________________________
		// 		get next batch number
		//________________________________________________________________________________
		
		t_batch_num batch_num (pg_db);
		if (args.use_batch)
		{
			if (!batch_num.get(args.job_category, args.batch))
				break;
		}
		else
		{
			if (!batch_num.get_next(args.job_category))
				break;
		}
		
		args.VERBOSE() << "batch=" << batch_num.text() << std::endl;
		
		
		pg_db.BeginTransaction();
		//________________________________________________________________________________
		// 		get gene pairs
		//________________________________________________________________________________
		t_gi_pairs 	gi_pairs;
		if (!batch_num.get_gene_pairs(gi_pairs, args.CERR()))
			return 1;
		args.VERBOSE() << "Retrieved gene pairs" << std::endl;
		
		//________________________________________________________________________________
		// 		get gene sequences, pis etc.
		//________________________________________________________________________________
		t_array_gi_seq array_gi_seq;
		t_map_gi_index map_gi_index;
		if (!get_gene_sequences(pg_db, batch_num.text(), array_gi_seq, map_gi_index, args))
			return 1;
		args.VERBOSE() << "Retrieved gene sequences" << std::endl;
		
		//________________________________________________________________________________
		// 		align all genes in pairs
		//________________________________________________________________________________
		deque<t_align_error>	alignment_errors;
		if (!align_all_pairs(pg_db, array_gi_seq, map_gi_index, gi_pairs,	
								blast_options_seg, blast_options,
								cnt_blast_comparisons, args, alignment_errors))
			return 1;
		
		if (!write_alignment_errors(pg_db, alignment_errors, args))
		{
			return 1;
		}
		
		//________________________________________________________________________________
		// 		batch succeed, free batch number
		//________________________________________________________________________________
		batch_num.commit();
		pg_db.Commit();
		
		++cnt_batches;
		if (args.use_batch)
			break;
	}
	}
	catch(...)
	{
		return 1;
	}

	args.VERBOSE() << "All batches finished\n";
	args.VERBOSE() << timer.pretty_reset() << " for " << cnt_batches << " batches\n";


	return 0;
}

